Quantum Dynamics of Three Coupled Atomic Bose-Einstein Condensates 



K. Nemotoll, C. A. Holmes*, G. J. Milburnt and W. J. Munro^ 

' Centre for Laser Science, Department of Physics, 

* Department of Mathematics, 
The University of Queensland, QLD 4072 Australia 
(February 1, 2008) 

The simplest model of three coupled Bose-Einstein Condensates (BEC) is investigated using a 
group theoretical method. The stationary solutions are determined using the SU(3) group under 
the mean field approximation. This semiclassical analysis using the system symmetries shows a 
transition in the dynamics of the system from self trapping to derealization at a critical value for 
the coupling between the condensates. The global dynamics are investigated by examination of the 
stable points and our analysis shows the structure of the stable points depends on the ratio of the 
condensate coupling to the particle-particle interaction, undergoes bifurcations as this ratio is varied. 
This semiclassical model is compared to a full quantum treatment, which also displays the dynamical 
transition. The quantum case has collapse and revival sequences superposed on the semiclassical 
dynamics reflecting the underlying discreteness of the spectrum. Non-zero circular current states 
are also demonstrated as one of the higher dimensional effects displayed in this system. 



I. INTRODUCTION 

The recent creation of neutral atom Bose-Einstein 
condensates (BEC) Jl|-^| stimulated theoretical research 
aimed at understanding this new state of matter. Models 
of two coupled BECs in a two mode approximation are 
considered a tractable system when total particle num- 
ber is conserved and eigenstates of the two well system 
are labeled by the particle number difference between the 
wells. Two coupled BECs in a symmetric double- well po- 
tential has been analyzed with the use of the SU(2) group 
|^],fr| to show the dynamical transition from self-trapping 
in one well to delocalized oscillation through both poten- 
tial wells due to the nonlinear particle interaction. This 
model in the weak coupling region has been further shown 
to demonstrate 7r-phase oscillations ||, while a semi- 
classical functional expression for the three-dimensional 
Josephson coupling energy have been derived ||. This 
model, however, can be considered as a special case be- 
cause of its simplicity and low-dimensionality. Any ex- 
tension of this model significantly increase the nonlinear- 
ity while higher dimensional effects increase the complex- 
ity of the model structure. 

These more complex systems are of interest as the 
richer dynamics and model structure allows us to treat 
quantum states with non-zero currents, for instance, fn 
the limit of large mode number n, Bose-Hubbard type ap- 
proaches are useful using a mean field approximation [|10| . 
However systems with intermediate numbers of modes, 
2 < n < 100, are complex and models must exploit sys- 
tem symmetries in order to obtain solutions. The sym- 
metries of these groups are the SU(n) group symmetries, 
fn this paper we analyze a three coupled BEC system 
using the operator algebra of SU(3). This could be re- 
alized as three spatially separated atomic Bose-Einstein 
condensates (BECs) as illustrated in the right corner of 



Fig. [I]. Alternatively it could describe a three conden- 
sates, occupying a single trap and distinguished by three 
internal hyperfine atomic states. The spatially separated 
system could represent a BEC confined in a three dimen- 
sional trapping potential with three harmonic minima in 
the x — y plane. Tunneling is possible between all three 
minima. This symmetric triple- well system represents 
the simplest two-dimensional generalization of the one 
dimensional double-well B which allows states of non- 
zero angular momentum. If the nonlinear interaction be- 
tween the atoms is not too large, see || (that is the total 
number of atoms N is not too big), we can describe this 
system using a minimum of three Bose modes for the 
quantum field; 
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V>(«>2/»*) = ^2[cj(t)uj(x, y) + c](t)u j (x,y)*], (1) 

where the uj(x,y) are an appropriate set of orthonor- 
mal single-particle mode functions for this potential, and 
the annihilation and creation operators Cj , cj satisfy the 
equal time commutation relations 

M] = <%. (2) 

As in the double- well case [pi, we can choose approx- 
imate single particle mode functions which are the lo- 
calized single particle ground states for each of the three 
wells. We further assume that these three lowest localized 
states are sufficiently well separated from higher energy 
states, and that interactions between particles do not 
change this basic property of the system configuration. 
Finally we assume that the total number of particles N 
are conserved. These assumptions allow us to treat this 
model system in a three mode approximation. Hence, the 
many-body Hamiltonian describing atomic BECs can 
be written in terms of the mode operators as 
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ft = u^c]c j +Q CjCk + X^cjctcjCj, (3) 

where o> is the mode frequency, x(< 0) is the two par- 
ticle interaction strength and f2(< 0) is the tunneling 
frequency. The condition, \ 5~ 0, corresponds to atoms 
with attractive interactions. This is the Hamiltonian of 
our model system in this paper. 



small number of atoms JlJ]. For attractive forces, the 
ground state is three-fold degenerate and in the occu- 
pation number representation (\m, n, N — (m + n)) = 
\m)i (g> \n)2 <S> \N — (m + these states are 

Id) = |0,0,7V), |ea) - |^,0,0), |e 8 ) = |0,iV,0). (7) 

For all these states the ground state energy is E = 
2 X N 2 /3. 



II. SU(3) GROUP APPROACH 

This section shows the group theoretical treatment of 
the system with the Hamiltonian specified in (||). In or- 
der to describe this system with SU(3) generators, we 
extend the Schwinger boson method [Ol and we define 
the eight generators of SU(3) {Zk,Yk, Xk} as 



c\ci - c\c 2 
X 2 = \ic\cx +4c2-2cjc 3 ) 



Y k = i(c\c : 



Z k = c\cj 



- c]c fe ) 



(4) 



where k = 1, 2, 3 and j = (fc + 1) mod 3. We note that 
the two operators Xi commute with each other. X\ and 
X 2 represent particle distributions projected on the y 
and x axes respectively in the right corner of Fig. |l|, and 
corresponds to angular momentum in this sys- 
tem. The most important relation which the generators 
satisfy is the Casimir invariant of SU(3), AN(N/3 + 1), 
where N is the total number operator, N = Y2j=i 



The operator algebra implies three further important 
identities; 



— + X2 + XAI— + X2-XA + — + 2x 2 

= y? + z\ 




= Y 2 2 + Zi (5) 
AN 




= y 3 2 + zi 

With the use of SU(3) generators, we can represent the 
Hamiltonian (^) in the form 

H = n(Z 1 + Z 2 + Z 3 ) + ^(X* + 3X$). (6) 

Here we ignore constant terms involving the conserved 
total number of particles N which do not change the dy- 
namics of the system. 

We now specifically consider the case with il — and 
X < 0. Such condensates are necessarily limited to a 



III. SEMICLASSICAL DYNAMICS 



We treat the three coupled BEC model using the semi- 
classical mean-field approximation. Ignoring correlations 
between all operators and taking expectation values con- 
verts the eight operator differential equations for the 
SU(3) generators in the quantum system into eight dif- 
ferential equations for real variables in the semiclassical 
system. It is, however, very difficult to analytically solve 
the full eight dimensional equations. This leads us to 
consider a specific simplified subspace of the set of eight 
equations by taking the symmetric condition, x\ = 0. 
The nature of this condition will be explained in subsec- 
tion A below as to see the dynamics clearly we first need 
to scale the semiclassical variables. 

The expectation values of generators are distinguished 
by their subscripts, while the expectation value of the 
total number operator is N . It is convenient to scale all 
the semiclassical averages by N. Thus we define, 
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(8) 



The equations of motion can be derived from the Hciscn- 
berg equations of motion of the Hamiltonian (Q), by fac- 
toring all higher order moments. 

The three degenerate ground states for £1 = can be 
associated with particular initial conditions in the semi- 
classical limit. To see this we first note that if we take 
matrix elements of both sides of the three operator iden- 
tities in Eq. (||) we find that 



(e*\Y? 



Zj\ 



(1 



Sij)2N. 



(9) 



Using the commutation relations for {Zi,Yi, Xi} and the 
corresponding uncertainty relations with respect to the 
ground states it is possible to show that the ground state 
variances in {Zi : Y{\ scale as N for N 1. In physical 
terms this means that the relative fluctuations in these 
variables goes to zero as N — > 00, as expected for a semi- 
classical limit. This indicates that in the semiclassical 
limit we may approximate 



{%) IN* = 4 

(Y*)/N 2 K(Ytf/N 2 = y* 



(10) 
(11) 
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The resulting semiclassical equations will be analyzed 
from two perspectives. We first show a dynamical transi- 
tion between self-trapping and derealization when initial 
conditions are given by Eq. (^) . This analysis determines 
the critical point for the dynamics transition. Secondly 
we consider the dependence of the stable points on the 
ratio of the coupling between condensates to the particle- 
particle interaction strength, and show bifurcations in the 
set of the stable points. 
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Using these relations in Eq. (||) we can construct semi- 
classical correspondences for each of the ground states. 
This is shown below 
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It is easy to verify that these curves are invariant under 
the semiclassical dynamics with O = 0. As each of the 
ground states is equivalent, up to a rotation in the phase 
space, we will now restrict the discussion to |ei) (a con- 
densate localized initially in the third well) without loss 
of generality, and examine the dynamics when CI ^ 0. 

Use of the initial state |ei) naturally restricts the dy- 
namics due to system symmetries, allowing us to study 
an analytically solvable subsystem. The initial state \e\) 
and the Hamiltonian (^) are symmetric to permutations 
of wells 1 and 2, then the resulting dynamics also satisfy 
this symmetry. This ensures that X\ = for all time 
in the semiclassical limit, which we previously referred 
to as the symmetric condition. This symmetric condi- 
tion specifies the subspace in which the reduced system 
dynamics lies, giving the additional conditions, yx = 0, 
?/2 = — 2/3, and z 2 = £3. The semiclassical dynamics is 
governed by the following four dimensional system 



x 2 = -2Cly 2 

i/2 = fi(3x 2 + z-y — z 2 ) - 3xNz 2 x 2 
Z\ = -2fh/ 2 

Z2 = Cly 2 + iX N X2V2- 



(12) 



The initial state \e\) gives the semiclassical system the 
initial conditions 



x 2 (0) = -2/3 
Vi(0) = 2((0) = 0. 



(13) 



The semiclassical dynamics governed by the equations of 
motion ( fl2"| ) is numerically shown in Fig. |l]. 



FIG. 1. Phase space orbits of the semiclassical dynamics 
projected on the 1/2 — £2 plane, for various values of r. The 



total number of atoms iV in the system was 5 x 10 



= 1/3 



is the critical value for localization. The transition from lo- 
calized to delocalized dynamics is apparent. The sub-figure in 
the right upper corner is a schematic representation of three 
spatially distinct Bose-Einstein condensates located at the 
minima of a potential with triangle symmetry. The tunneling 
coupling constant between all the wells is equal and nonzero. 

Fig. [I] shows a transition from localization to global 
oscillation, and this dynamics change can be explained 
by the stable point analysis. This reduced system is inte- 
grable as these equations satisfy two constants of motion; 



Q,(zi + 2z 2 ) + 



3 X Nx 2 2 



3x 2 2 + 2(y| + z 2 2 ) + z\ 



H/N 
4/3, 



(14) 
(15) 



which correspond to energy and total number conserva- 
tion respectively. The latter constant ( |T5| ) follows from 
two stricter constraints 



X2 



2 [x 2 + - 



2,2 1 

V2 + z 2 = 2> 



(16) 



(17) 



which are derived from Eq. (|^) with the symmetric con- 
dition in the semiclassical limit. 

The equations of motion with the chosen initial con- 
ditions may be solved explicitly as a function of X2 and 
depend on only one parameter; the ratio of the coupling 
constant $1 to the interaction strength \, 



r=^(> 0). 



(18) 



The solutions are, 



Zi(t) 

*2(t) = - 



x 2 (t) + 
x 2 (t) 



3x 2 (t) 2 
4r 



1 
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(19) 



V2{tf = f(x 2 (t)), 
where 
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fMt)) = -Ur)-(x 3 (t) + ^Y^ 

+y (» 2 (t) + 0(2-i)+B(r), (20) 

and where the integration constants are given by 
Ar (r-2) 2 



a{r) 
B{r) 



3 12r 

2 

*(2r-l)(r-2). 



(21) 
(22) 



The solutions are oscillatory and only exist for —2/3 < 
x 2 < 1/3, and are strongly dependent on the roots of the 
function fix), f(x 2 ) is fourth order in terms of x 2 and 
can have up to four real roots. For increasing r, the root 
structure of f{x%) changes. To really understand these 
changes it is instructive to look at the fixed point struc- 
ture of system (12). If all the derivatives of the equations 



of motion (|12j) are zero then, 
V2 = 0, 



-3x 2 — Zi+ z 2 + -z 2 x 2 = 0. 
r 



(23) 



Taking the positive root for Z\ in the constraint (17) 
z\ = x 2 + 2/3, the second equation can be written 



3 2 
4x 2 + z 2 + -z 2 x 2 - t: = 0. 
r 3 



(24) 



So the fixed points are the crossing points of this equation 
with the constraint, 
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6/ ' " z 2' 



2 ( x 2 + - ) + zf 



(25) 



obtained by applying y 2 = to constraint ( |l7|) . Together 
these result in the following fourth order ploynomial in 

z 2 
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To I z 2 - 77 



2 Ar 1 
-(1 - 4r)z 2 - 2r(l - r)z 2 H 1=0 



Fig. |^ shows the curves (^4|) and (|2^) in the x 2 -z 2 plane 
for various values of r. The first factor of ( |26| ) gives us 
the point A (x = 0, z 2 = §)■ For larger r only one branch 



of the hyperbola (|24|) intersects the constraint and there 
are just two elliptic fixed points, one at A and one in 
the third quadrant. The solutions with initial condition 
x 2 = — | are far away from the elliptic fixed points and 
perform large delocalised oscillations on the sphere (p~7h . 
A projection of one such solution is shown in Fig. RTfor 
r = 0.52. As r is decreased a second branch touches the 
constraint at C Fig. ||. 




FIG. 2. The r dependence of the stable points. Each hy- 
perbolic curve shows equation (^i|) for different r, while the 
ellipse represents constraint (|25|). Point A is a stable point for 
any r. Point B shows when the equation collapses into 
the two lines, and point C is the unique point of tangential 
contact to the ellipse. 

This occurs when the second factor in (^6|) has a dou- 
ble root at r = 0.507425. Two new fixed points, a saddle 
and a center, then appear and move apart as r is de- 
creased away from 0.507425. The stable and unstable 
manifolds of the saddle intersect, forming a lopsided fig- 
ure eight like separatrix with the new elliptic fixed point 
in the one lobe and the elliptic fixed point at A in the 
other. When r is | the hyperbola reduces to the lines 
z 2 = \,x 2 = — i and the solutions in (x 2 ,y 2 ) space are 
symmetrical about x 2 = — ^. If r is further decreased 
the separatrix approaches our initial condition (|l3|). At 
r = | the initial condition actually lies on the separatrix. 
This occurs when the unstable fixed point lies on the so- 
lution curve (pit) , which amounts to looking for a double 



root of f(x 2 ). If we let u — x 2 



4r 



then 



2u 



+ —(2r-l) + B(r). 
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and we require / = and /' = 0. /' = gives, 
u 2 = y (a{r) - v/4a(r) 2 - 35(r)) , 



(27) 



(28) 



■(wMijbh when substituted into / = gives only one real 
solution at r* = 1/3, x 2 — 0. Now for r > | the solu- 
tions lie within the separatrix and so as r is decreased 
through -| the oscillations suddenly reduce to half their 
former size as can be clearly seen in Fig. ^. 

In the physical space of the potential then, the con- 
densate remains localized at the bottom of the first 
well,a;2 = — §, for r = . As we increase r from zero, 
x 2 (t) begins to oscillate. In the x 2 — y 2 plane, the point 
is replaced by small oscillations, within the separatrix, 
which pass through x 2 = —2/3. This is referred to as dy- 
namical localization. Even though the coupling between 
wells is present, the nonlinear interaction between parti- 
cles prevents the condensate from moving away from its 
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initially localized state. Eventually for the critical value 
of r= g the orbit in the X2 — J/2 plane extends across both 
half planes for positive and negative X2 values and the 
condensate is no longer localized. 

In Fig. [l] we show the phase space orbits of the semi- 
classical dynamics projected onto the X2 — IS2 plane for 
various values of the parameter r. The initial conditions 
correspond to a condensate localized in well 3. The tran- 
sition from localized to delocalized motion is seen when 
the orbit in the phase space extends into the X2 > half 
of the plane. In Fig. || we plot x% it) as a function of time 
for the condensate initially localized in well 3 with the 
same initial conditions as in Fig. |l| for two values of r, 
one above and below the critical value r* for localization. 
The transition from localized to delocalized oscillations 
is apparent. 
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3. Time evolution of X2 (t) in the semiclassical regime 
versus time below (Figure a) and above (Figure b) the crit- 
ical value r* = 1/3 for localization. Figure a) corresponds 
to the below threshold regime with r = 0.283 while figure b) 
corresponds to the above threshold regime with r = 0.506. 
The number of atoms is fixed at N = 50, and the time t is 
normalized by Q. 



quantum mean values decay due to the intrinsic quan- 
tum fluctuations in the number of atoms in each individ- 
ual well, while the total particle number in the system 
remains fixed. The collapses and revivals of the oscil- 
lations in the quantum system arise from the discrete 
nature of the eigenvalue spectrum for finite atom num- 
ber. Such phenomena were also observed in two coupled 
condensates [pi. 
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FIG. 4. The quantum dynamics of (X2) versus time below 
(Figure a) and above (Figure b) the critical value r* = 1/3 
for localization. Figure a) corresponds to the below thresh- 
old regime with r = 0.283 while figure b) corresponds to the 
above threshold regime with r = 0.506. The number of atoms 
is fixed at N = 50, and the time is normalized by Q. 

One of the interesting higher dimensional features of 
this system is the existence of non-zero circular current 
states. For r — 0, the system has the three-fold degener- 
ate ground states, |ei), |e2), and |e3), and superpositions 
of those states create another set of orthonormal ground 
states given as 



IV. QUANTUM DYNAMICS 

In this section we treat the three coupled BEC model 
fully quantum mechanically and numerically calculate 
the time evolution of the particle distribution with the 
initial condition |ei). For fixed total atom number N, a 
suitable basis of the system Hilbert space is the number 
eigenstates |n, m,N — (n + m)) which are simultaneous 
eigenstates of the generators, X\ and X 2 of the Abelian 
sub-algebra of SU(3). The unitary evolution operator is 
given by U(t) = exp [— iHt] where H is specified by (||). 
The unitary evolution matrix is then computed and the 
initial state |^>(0)) = \e\) is then propagated forward in 
time. At each time step we compute the averages of (X\) 
and (A2). These averages show the particle distribution 
projected on the y and x axes in Fig. |l| respectively. How- 
ever because of the initial conditions for the state | -0(0)} 
the average of (X%) does not change and in fact remains 
zero. It is not considered further here. 

Fig. H shows the evolution of a condensate initially 
localized in state |ei) with number of atoms fixed at 
N = 50. For short times we see the same oscillations 
as in the corresponding semiclassical case, both below 
and above the critical value r». The oscillations of the 



l.9i> = -7= (|ei) + |e 2 ) + |e 3 » 

\92) = ^(e-^ /3 \e 1 ) + \e2)+e^ 3 \e 3 )) 

\9s) = j= (e^ /3 |ei) + e^ 3 \e 2 ) + |e 3 >) . 



(29) 



These states are invariant under 2tt/3 rotations due to 
system symmetries, as discussed in UM. Taking the state 
1 92), we examine the quantum dynamics of the angular 
momentum Y s = Y\ + Y 2 + Y 3 . (The state \g 3 ) gives 
the same dynamics, however since these states are anti- 
symmetric to each other, clockwise motion in | (72) cor- 
responds anti-clockwise in 1 53 } - ) For r = 0, the state 
1 52) is the ground state, and the average angular mo- 
mentum remains zero, though for non-zero r, a non-zero 
average angular momentum appears as seen in Fig. |^. 
This shows quantum dynamics of the average angular 
momentum normalized by the total number iV for two 
different values of r. For small r the non-zero angular 
momentum does not develop into any stable circular mo- 
tion (Fig. ||a), while circular motion can be established 
for large r (Fig. |^b), and becomes increasingly stable for 
larger r. 
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FIG. 5. The quantum dynamics of (Y a )/N versus time with 
the initial state |<?2), where Y s = Yi + Yi + Y3 is the angular 
momentum. The angular momentum for small r = 0.5 (Fig- 
ure a) fluctuates, while non-zero circular motion appears in 
Figure b) for larger r — 1.3. The number of atoms is fixed at 
N — 50, and the time is normalized by Q. 



V. DISCUSSION AND CONCLUSION 

In this paper we have shown how the generators of 
SU(3) can be used to describe the quantum and semi- 
classical dynamics of three symmetrically coupled atomic 
Bose-Einstein condensates. By taking expectation val- 
ues of the Heisenberg equations of motion and factoring 
all higher order moments we can derive the semiclassi- 
cal mean field equations. The nonlinear terms arising 
from hard collisions lead to a dynamical bifurcation in 
the semiclassical dynamics as the tunneling strength is 
increased, reflecting a transition from localized dynamics 
to tunneling currents. The r dependence of the system 
is much more complicated than that found for two cou- 
pled BECs. The r dependence of the stationary solutions 
in this paper constitutes only a small part of the total 
complexity of this system and only some of the higher 
dimensional effects. Our quantum treatment verified the 
dynamical transition found in the semiclassical analysis. 
This three coupled BEC model is the simplest model to 
have non-zero circular current states, and we have shown 
non-zero circular motion can appear given appropriate 
initial conditions. 

The analysis of this paper focused on comparing the 
semiclassical dynamics with the evolution of quantum 
mean values. This restricted analysis naturally suggests 
further examination of the dynamics of full quantum 
states. However in order to examine full quantum states, 
it is necessary to use more powerful group theoretic tools 
and this will be the subject of a future paper. 
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